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Similar evolutionary variational inequalities appear as convenient formulations for 
continuous quasistationary models for sandpile growth, formation of a network of 
lakes and rivers, magnetization of type-II superconductors, and elastoplastic defor- 
mations. We outline the main steps of such models derivation and try to clarify 
the origin of this similarity. New dual variational formulations, analogous to mixed 
variational inequalities in plasticity, are derived for sandpiles and superconductors. 



O 



> 



cd 



Key words: critical states, variational inequalities, sandpiles, superconductors 



(N) ' PACS: 02.30.W 

O 
O 



Introduction 
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C . 

Spatially extended dissipative systems have recently attracted much interest 

among physicists. These systems may have infinitely many metastable states 

but, driven by the external forces, often organize themselves into a marginally 

^ stable "critical state" and are then able to demonstrate almost instantaneous 

long-range interactions between their distantly separated parts. 

Modifications of a crude cellular automata model of sandpile [1] have been used 
by many authors for simulating such systems behavior and it was sometimes 
doubted whether the models based on differential equations can in principle 
be employed: the relaxation in continuous models is expected to be a smooth 
process evolving in time, while, e.g., cellular automata models are able to 
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mimic sand avalanches as sudden catastrophic events. Nevertheless, continu- 
ous models allowing for long-range interactions, hysteresis, metastability, and 
avalanches have been derived for sandpiles [2,3,4], river networks [2], type-II 
superconductors [5,6,7]. Although these are dissipative systems of a differ- 
ent nature, their continuous models are equivalent to very similar variational 
or quasivariational evolutionary inequalities, the formulations convenient for 
both the numerical simulation [8] and theoretical study [9]. Much earlier, sim- 
ilar variational formulations have been derived for various plasticity problems 
(see [9]). 

The aim of this work is to outline the main steps of such models derivation 
and to clarify the origin of this similarity. We present the simplest version 
of each model and try to avoid mathematical details but give references to 
rigorous proofs, numerical simulations, and possible extensions of the consid- 
ered models. We also derive new dual variational formulations in terms of the 
conjugate variables for sandpiles and superconductors. The dual problems are 
similar to mixed variational inequalities in plasticity [10]. Well-posedness of 
these new problems is yet to be investigated. 



1 Variational formulations 

1.1 Sandpiles 

Let a cohesionless granular material, characterized by its angle of repose a, 
be poured out onto a rough rigid surface y = ho(x), where y is vertical and 
x G VL C M 2 . We find the shape of a growing pile, y = h(x, t). 

Assuming the flow of granular material down the slope of the pile is confined to 
a thin boundary layer and the bulk density of material in the pile is constant, 
we can write the mass conservation law in the form dth + V ■ q = w, where q is 
the horizontal projection of the material flux and w(x,t) > is the intensity 
of the source of material being poured onto the pile. Neglecting inertia, we 
suppose that the surface flow is directed towards the steepest descent, q = 
—mWh, where 

m(x,t)>0 (1) 

is an unknown scalar function. The conservation law assumes the form 

dth - V • (fflVA) = w. (2) 

The free surface initially coincides with the support, 

h(x,0) = ho(x), (3) 



and cannot lie below it, 

h(x,t)>h (x). (4) 

Wherever the granular material covers the support, the surface slope cannot 
exceed the material repose angle a, 

h(x,t)>h (x) — ► \Vh(x,t)\ <7, (5) 

where 7 = tan(o;). No surface flow occurs over the parts of the pile surface 
inclined less than at the angle of repose: 

\Vh(x,t)\ <7 — ► m(x,t)=0. (6) 

We assume for simplicity that there is a vertical wall at the boundary T of 
domain Q, hence 

md n h = on T. (7) 

The formulated model for pile growth contains two unknowns, the free surface 
h and an auxiliary function m, and it is difficult to deal with the equations 
and inequalities (l)-(7) directly. Fortunately, a more convenient variational 
formulation can be derived (here we follow [2], see [11] for a rigorous proof). 
Let us define, for every continuous function ip, a nonlinear operator 



B^) = l -(\V V \ 2 -M{i,)) 
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where 

{7 2 ifib(x.t) > h (x), 

max(7 2 , \Vh (x)\ 2 ) if ip(x, t) < h (x). 
We define also a family of closed convex sets l 

Kty) = {<p(x,t) \B4<p)<0}, 

denote by (u, v) the scalar product of two functions, and consider an evolu- 
tionary quasivariational inequality written for the pile surface alone: 

Find h G K(h) such that (d t h — w,tp — h) > for any (p G K(h), 

(8) 

h(x, 0) = ho(x). 



Theorem. A function h(x,t) is a solution of the quasivariational inequality 
(8) if and only if there exists m(x, t) such that the pair {h, m} satisfies a weak 
form of the problem (l)-(7). 



1 More exactly, K(ip) = {y? G L°°(0,T; W 1>00 (n)) | Bj,(ip) < a.e.}, see [11]. 



Proof. We formally rewrite the inequality (8) as an implicit optimization 
problem 

J h (h) = min J h (ip) 

(9) 

y? e K(h) 

where Jh(<p) = (d t h—w,ip) is a linear functional which depends on the solution 
h. Let us fix the function h in J^ and K(h) and derive the necessary and suffi- 
cient condition of optimality for (9) using the Lagrange multipliers technique 
([12], ch. 3, th. 5.1). Substituting then the function h into this condition, we 
obtain a similar condition for the problem with an implicit constraint: h is 
a solution of the quasivariational inequality (8) if and only if there exists a 
Lagrange multiplier m(x,t) > such that the pair {h, m} is a saddle point of 
Lagrangian, i.e., 

J h (h) + (m', B h (h)) < J h (h) + (m, B h (h)) < J h (ti) + (m, B h (ti)) (10) 

for arbitrary h! and m' > 0. The condition of supplementary slackness, 

(m,B h (h))=0, (11) 

is thereby fulfilled. 2 

Let h be a solution of (8). As follows from (10), the functional (d t h — w, h') + 
|(m, |V/i'| 2 — M(h)) has a minimum at the point h' = h. Hence, 

(d t h-w,x) + (m,Vh-Vx)=0 (12) 

for any test function x- This is a weak formulation of equation (2) with the 
boundary condition (7). Since h G K(h), condition (5) is satisfied, (6) follows 
from (11), and to show that {h, m} satisfies all model relations it remains only 
to check that h > h . Choosing 

h + (h - h)+ for < t < t , 

h otherwise, 

where z + means max(z, 0) and taking into account that ip G K(h), w>0we 
obtain 

1 /■ ,2 



0<(dth-w,il>-h)<--J {[ho(x) - h(x,to)] + y d&l, 
so the inequality (4) is proved. 

Let now {h, m\ be a solution to (l)-(7). By (4), (5) we have \Vh\ < 7 wherever 
h > ho and h = ho otherwise, hence h G K(h). To prove that h solves the 



2 As is explained in [11], to satisfy a constraint qualification hypothesis ([12], ch. 
3, (5.24)) we need to define B h : L°°(0,T; W 1 ' 00 ^)) -> C = L°°((0,T) x n). Hence 
m belongs to the dual space C and is a nonnegative Radon measure. 



quasivariational inequality, it is sufficient to show the inequalities (10) hold. 
It is easy to see that (11) is true, so the left of the inequalities (10) is fulfilled. 
Using the weak form (12) of equations (2), (7) we obtain, for x — h' — h, 

J h (ti) + (m,B h (ti)) - J h (h) - (m,B h (h)) = l -{m, \V{h' - h}\ 2 ) > 0. 

Thus the second inequality in (10) holds too, which completes the proof. 

We note that the auxiliary unknown, m, introduced into the pile growth model 
to fix the possible sand flux direction, turns out to be a Lagrange multiplier 
related to the equilibrium constraint upon the pile surface incline and is elim- 
inated in transition to the variational formulation. The multiplier depends in 
a non-local way on the surface and source and that is why instantaneous long 
range interactions over the critically inclined parts of the surface are possible 
in this model. Such a situation is typical also for other dissipative systems 
where the relaxation is fast and the assumption that all the dynamics occur 
at the border of stability is justified. 

If the support has no steep slopes, \Vh \ < 7, the set of admissible functions 
K becomes fixed (does not depend on the solution) and the inequality (8) 
becomes simply variational; in this case the existence and uniqueness of a 
solution have been proved [11]. It remains an interesting open problem to 
prove existence and uniqueness in the general quasivariational case. 

The variational formulation obtained is very convenient for numerical simu- 
lation of pile growth, see [13]. There are analytical solutions [2] exactly de- 
scribing the pile shapes in experiments [14]. Mathematically, the avalanches 
upon pile surface correspond to solutions with the jumps caused by sudden 
variations of the admissible functions set K due to local fluctuations of the 
repose angle, see [2]. Such discontinuous solutions of the variational inequal- 
ity have been studied in [4]. It has been also shown [15] that the mesoscopic 
BCRE model [16] for sand surface dynamics converges in the long scale limit 
to the inequality (8). In a continuous limit, stochastic cellular automata mod- 
els of sandpiles converge to a similar variational inequality with the anisotropy 
inherited from the cellular structure of these models [17]. 



1.2 Lakes and rivers 



Let now /i be the land surface, w the intensity of precipitation. We assume 
for simplicity that the water neither evaporates nor penetrates the soil but 
just flows down the slopes and accumulates into lakes at local depressions of 
the relief. The level of a lake rises until it reaches the divide of two basins. 
Then a river running out of the lake appears and transfers all additional water 



to another lake below. 

To model the evolution of this system of lakes and rivers, let us note that the 
balance equation 

d t h + V • q = w, 

in which q is the horizontal projection of water flux, remains valid. The free 
boundary h in this problem either coincides with h or, where it is higher, is 
the horizontal surface of a lake: 

h(x,t) > ho(x), h(x,t) > ho(x) — > V/i(x, t) = 0. 

Over the hill slopes, where h = ho, we assume again that the flux is in the 
steepest descent direction, 

h(x,t) = ho(x) — > q = —rnVh, 

where m(x, t) > is unknown. However, this is not true for the lakes, where 
h > h and V/i = 0. In fact, although the lake hydrodynamics are not trivial, 
the flow in the lake does not affect the free surface, and it can be shown [2] 
that the model relations above lead to the quasivariational inequality (8). 

Indeed, let 7 = and the pair {h, m} satisfies these relations. Then h G K(h) 
and for any function ip G K{h) we obtain 

(d t h-w,ip-h) = -(V -q,(p-h) = - i q n (ip - h) + / q-V(ip-h). 

Jo Jr Jo Jn 

The first integral on the right hand side is zero due to the boundary conditions. 
Gradients Vh and Wip are both zero wherever h > /i . Outside this domain 
h = h , q = —rnVho, m > 0, and |Vy?| < |V/io|- Therefore, 

q ■ V(y? - h) = -m(X7h ■ Vy? - |V/i | 2 ) > 0, 

the quasivariational inequality holds and determines the free surface evolution. 

This is, however, only a part of the solution needed: it is the water flux q = 
\q\, or, equivalently, the auxiliary variable m, which is of interest in most 
geomorphological and hydrological applications. Provided the free surface h 
is found, the water flux in the coincidence set h = ho can be determined, at 
least in some simple cases [2] , from the water balance equation which takes in 
this set the form 

V7 ( V M 

V |v/io|y 

Generally, this is a difficult task and a different approach to water flux calcu- 
lation is desirable. Below, we consider an alternative approach to determining 
the conjugate variables for variational inequalities. 



Note that if 7 = the quasivariational inequality (8) is no more equivalent to 
the sandpile model (l)-(7) in which q — if Vh = 0. This equivalency breaks 
down because, if 7 = 0, the constraint qualification hypothesis [12] is not true: 
there exists no ip such that B h (ip) < 0. 



1.3 Superconductors 



Phenomenologically, the magnetic field penetration into type-II superconduc- 
tors can be understood as a nonlinear eddy current problem. In accordance 
with the Faraday law of electromagnetic induction, the eddy currents in a 
conductor are driven by the electric fields induced by time variations of the 
magnetic flux. Let the superconductor occupy a domain QcK 3 and to = M?\Q 
be the outer space. We denote by T the common boundary of these domains 
and assume n, the unit normal to I, is directed outside Q. 

Omitting the displacement current in Maxwell equations and assuming the 
magnetic permeability of the superconductor is equal to that of vacuum and 
scaled to be unity, we obtain the following eddy current model, 

d t h + V xe = 0, 

xeM 3 , t > 0, (13) 

V x h = j + j e , 

with the initial condition h\ t =o = ho(x) such that V • h$ = 0. Here j is the 
current in the superconductor (j = in 00), and j e is the given external current 
having a bounded support supp{j} C 00 and satisfying Vj e = 0. Additionally, 
in the conductive domain Q a current-voltage law has to be postulated. 

In an ordinary conductor, the vectors of the electric field and current density 
are related by the linear Ohm law. Type-II superconductors are instead char- 
acterized in the Bean critical-state model [18] by a highly nonlinear current- 
voltage relation which gives rise to a free boundary problem. The problem has 
been solved mainly under the assumption that the electric field and current 
density are parallel (see, e.g., [7] and the references therein). Then e = pj, 
where the a priory unknown effective resistivity p(x, t) > characterizes the 
energy losses accompanying movement of magnetic vortices in a superconduc- 
tor. It is assumed in the Bean model that the current density j cannot exceed 
some critical value j c and, until this value is reached, the vortices are pinned 
and the current is purely superconductive: 

\j(x,t)\<j c , \j(x,t)\<j c --+e(x,t) = 0. (14) 



The simplest geometric configuration is that of a long superconductive cylin- 
der having a simply connected cross-section Q C M 2 and placed into a non- 



stationary parallel uniform external magnetic field h e (t). In this case the most 
convenient variational formulation can be derived for the magnetic field in the 
superconductor. The current density j induced by the external field variations 
is parallel to the cross-section plane and produces the magnetic field hi(x,t), 
parallel to h e and equal to zero on T. The problem is two-dimensional. De- 
noting by hi and h e parallel to the cylinder axis components of hi and h e , 
correspondingly, and using the standard notations curlv = d Xl V2 — d X2 Vi, 
curlv = (d X2 v, —d xi v), we rewrite the model (13) as 

d t (hi + h e ) + curl e = 0, curl hi = j, x G f2, t > 0. (15) 

Since |V/ij| = \curlh{\ = \j\ < j c , hi(x,t) should belong to the set 

K = {<p(x,t) | \Vcp\<j c , ^|r = 0}. (16) 

Multiplying the first of equations (15) by ip — hi, •p G K, integrating, and 
using the Green formula we obtain (d t {hi + h e }, tp — hi) = (e, j) — (e, curl (p). 
Taking the Bean current-voltage relations into account we get 

(e,j) = (|e|,|j|) = (\e\,j c ) > (\e\,\curlip\) > (e,curlip) 

and arrive at the variational inequality for the induced field hf 



Find hi G K such that (d t {hi + h e }, ip — hi) > for any ip G K, 
h i (x,0) = ho(x)-h e (0). 



(17) 



In the general three-dimensional case the /i-formulation of the Bean model can 
also be derived [6] but then the magnetic field should be determined in the 
whole space. A similar variational formulation in terms of the current density 
we derive now is probably more convenient. 

Provided that no current is fed into a superconductor by electric contacts, 
V • j — in n, j n = on T. Let us define the set of admissible current 
densities in Q, 



j eK = { cp(x, t) 



V ■ cp = 0, \<p\ < j c in Q 
cp n = on T 

express the electric field via the vector and scalar magnetic potentials [26], 

e + d t A + V^ = 0, 

and exclude the scalar potential by multiplying this equation by <p — j and 
integrating: (e + d t A, <p — j) = for any cp G K. Just as above, if e and j are 
parallel and satisfy the Bean model relations (14), the inequality (e, (f—j) < 
holds for every admissible test function cp. Hence, (d t A, ip — j) > 0. 



Up to the gradient of a scalar function, determined by the gauge and elim- 
inated by scalar product with the test functions, the vector potential is a 
convolution of the Green function of Laplace equation, G = jy-r, and the 
total current: 

A = G*{j+j e }. 

We arrived at the evolutionary variational inequality with an "implicit" deriva- 
tive in time: 

Find j G K such that 
(G * d t {j + j e }, <p-j)>0 for any cp e K, (18) 

j\t=o = jo(x), 

where jo — V x ho\u G K is a given initial current density distribution. 

Experiments on hard superconductors are often performed on thin flat sam- 
ples, and we present also a scalar version of this variational inequality for thin 
films in a perpendicular uniform external magnetic field. We now assume that 
it is the sheet current density, obtained by integration of bulk current density 
across the film thickness and also denoted j(x,t), x G Q C 1R 2 , obeys the 
Bean's current-voltage relations. This current density should also satisfy the 
conditions 

divj = 0mn, j' n = on r, (19) 

where div is the two-dimensional divergence. For simplicity, we assume the 
domain Q is simply connected. Due to conditions (19) there exists a stream 
function h(x, t) such that curl h = j in Q and h — on I\ Although h is 
not the induced magnetic field as in the case of a long cylinder in a parallel 
field, since \curl h\ = \Vh\ this function belongs to the same set K of admis- 
sible functions (16). Let j' be another vector function satisfying (19) and the 
condition \j'\ < j c in Q, and tp G K be the corresponding stream function. 
The external vector potential A e , corresponding to the uniform perpendicular 
magnetic field h e (t), can be chosen parallel to the film; then curl A e = h e (t). 
Substituting into the inequality (18) the curls of h and <p instead of the cur- 
rent and test function, correspondingly, using the Green theorem, and taking 
into account that curlu ■ curlv = Vw • Vi> we obtain a scalar variational 
inequality in terms of the stream function: 

Find h G K such that 
a(d t h, (p — h) + (d t h e , p — h) > for any p G K, (20) 

h(x, 0) = h (x), 

where a(u,v) = J n J n Vu(x) ■ Vv(x')/{4ir\x — x'\} dxdx' . 

The existence and uniqueness of solutions to (18) and (20) were proved in [6] 
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and [19]; see [20] for the numerical solution of (20). It has been shown in [6] that 
the effective resistivity p, excluded in transition to the variational formulation 
(17), is a Lagrange multiplier related to the current density constraint. Sim- 
ilar variational formulations may be derived for much more general current- 
voltage relations (see, e.g., [5,7,21]) and present a very convenient description 
of hysteretic magnetization typical of hard superconductors. In particular, the 
critical current density j c depends usually on the magnetic field [23]. Then 
the set of admissible functions K depends on the unknown solution and the 
inequalities become quasivariational [6]. The power law |e| = e (\j\/j c ) p is 
often employed instead of the Bean's current voltage relation to account for 
the creep of magnetic flux [22]; as p — » oo, such model converges to the Bean 
model [19,24]. Thermal fluctuations in a superconductor may cause avalanches 
of magnetic vortices resembling sand avalanches [25] ; in the Bean model these 
avalanches correspond to discontinuous solutions of variational inequality (18) 
with the jumps induced by local fluctuations of the critical current density. 



1.4 Elastoplastic solids 



The variational inequality formulation for models in perfect elastoplasticity is 
well known [9]. We briefly present this formulation to underline its similarity 
to the variational formulations above. 

Let an elastoplastic body occupy the domain Q C 1R 3 and the conditions of 
equilibrium, 

fg+ff = 0, fxxg+fxxf = 0, 
Ja Jr Jn Jr 

hold for the given body force g and surface traction /. The stress tensor cr 
should satisfy the local equilibrium conditions (the usual summation conven- 
tion is implied) 

Vijj + 9i = in Q, (Jijiij = fo on T. (21) 

Under the assumption of small strain we have 

e v = 2^' + "**)' ( 22 ) 

where u is the displacement vector field, Uij = dui/dxj. 

It is assumed that the strain tensor can be presented as a sum of elastic and 
plastic components, e=e+p, where the elastic component obeys the linear 
Hooke's law, e„ = A^m <Jki- The plastic part is governed by an incremental flow 
rule, p = jQpdt, in which p = d t p is determined as follows. For a prescribed 
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convex yield function J- {a), it is postulated that the stress tensor everywhere 
satisfies 

F(&) < (23) 

and that 

where X(x, t) > is the deformation rate such that T((r(x,t)) < — ► 
X(x,t) = 0. Let the admissible set K be the set of stress tensors satisfy- 
ing (21) and (23) and r G K. Using the strain-displacement relation (22), the 
Green formula, and the equilibrium conditions (21), one can show [27] that 
(tij,Tij -&ij) =0, hence 

(Aijki &kl, Tij - Gij) + (X d? '/d(Jij, Tij - (Tij) — 0. 

The second term here is nonpositive. Indeed, if J-(cr(x, t)) < then X(x,t) = 
0. Otherwise jF(cr(x,t)) = and, because T is convex, ^({l — 9}a(x,t) + 
9r(x,t)) < for any 9 e [0, 1]. Therefore, 



&W-9W 



e=+o vcTij 



{Tij - (Tij) < 0. 



We arrived at the variational inequality: 

Find a G K such that (A^ki &ki, r^ — a^) > for any r G K, 

<t(x, 0) = (Tq(x). 



2 Dual formulations for conjugate variables 



Although a common feature of the considered mathematical models is the 
presence of conjugate variables, the variational inequalities above are writ- 
ten for only one of them: the free surface in the model for sandpile growth, 
current density in critical-state superconductivity problems, stress tensor in 
elastoplasticity problems. The dual variables, i.e, the surface flux, the electric 
field, and the strain, correspondingly, have been eliminated in transition to the 
variational inequalities. These inequalities can be solved efficiently, however, 
knowledge of the primary variables generally does not make determining the 
dual ones easy because the constitutive relations are multivalued. 

In elastoplasticity with hardening, a dual variational formulation for strain has 
been derived and comprehensively studied by Han and Reddy [10]. Mathemat- 
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ically, this problem takes the form of the so-called mixed variational inequality, 

Find u : [0, T] — ► V such that for almost all t 
a(u{t),v - d t u(t)) - (f(t),v - d t u(t)) + <f>[y) - <f>(d t u{t)) > for all v e V, 

u(0) = u , 

where V is a Banach space, a(., .) is a symmetric bilinear form, f(t) G V is 
a linear functional, and 0(.) is a convex, positively homogeneous, nonnegative 
functional on V. This problem is not a standard variational inequality although 
it resembles the parabolic variational inequalities of the second type [8] . 

Strain formulation for the perfect plasticity problem considered above turns 
out to be much more complicated because, without hardening, the arising 
problem is not coercive in the usual Sobolev spaces (coercive in a non-reflexive 
Banach space) and the solution has to be sought in the space BD(Q) of func- 
tions of bounded deformation ([28], see [29] for a review of recent results). 
Physically, the difference is manifested in the ability of perfectly plastic mate- 
rials to form slip surfaces on which the tangential component of displacement 
is discontinuous. 

Similar difficulties arise in dual variational formulations of other critical-state 
problems described in the first part of our work and this, indeed, seems to be 
physically meaningful. Thus, the continuity of only the normal component of 
surface sand flux follows from the mass conservation law in the pile growth 
model. Also, according to Maxwell equations, only the tangential component 
of the electric field has to be continuous. Mathematically, the corresponding 
problems are not coercive in reflexive Banach spaces. 

Below, we derive formally variational formulations in terms of the dual vari- 
ables for two critical-state problems where the primary formulation is a varia- 
tional (not quasivariational) inequality. The questions of existence, uniqueness, 
and numerical approximation of these problems need further investigation; we 
are going to consider these questions in a separate publication [30] . 



2. 1 Surface flux in the model of pile growth 



Determining the flux of granular material pouring down the free surface of a 
growing pile is necessary, e.g., if the material is polydisperse and it is needed 
to predict the resulting distribution of different species inside the pile (see 
[13]). Let us assume the initial support h in the pile growth model has no 
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steep slopes, |V/io| < k in fi. In this case the model (l)-(7) can be written as 

d t h + V-q = f, 

y J ' (24) 

h\ t =o = h , q n \r = 0, 

where the flux q has the direction of —Vh and the following flux-slope relation 
holds: 

\Vh\<k, \Vh\<k->q = Q. (25) 

As it was shown above, the free surface h(x, t) can be sought as a solution of 
an evolutionary variational inequality. 

To derive a variational formulation of this model in terms of the surface flux, 
let us define 

V = H (div; O) = {cp G L 2 (n) \V-<pe L 2 (n), <p n \ r = 0}, 

assume that q and h satisfy the model relations (24)- (25), and choose an 
arbitrary test flux q £ V. Using the constitutive relations (25) we obtain 

Vh-(q-q)> -\Vh\\q\ -Vh-q= -\Vh\\q\ + k\q\ > -k\q\ + k\q\. 

Hence, 

{Vh,q-q)><i>{q)-<j>{q), 

where (p(q) = k J n \q\. Since (Vh,q — q) — — (h, V ■ {q — q}), we have 

(f>(q)-cf ) (q)-(h,W-{q-q})>0. 

Let us define u = J f qdt. Then d t u = q and, from (24), V • u = —h + /i + 
/o / dt. We finally arrive at the following variational problem: 

Find u : [0, T] — > V such that for any q G V and almost all t 
(V • u, V • {q - d t u}) - {T, V • {q - d t u}) + <f>{q) - <f>{d t u) > 0, (26) 

and u\ t =o = 0, 

where T = ho + Jq f dt. Since the problem is not coercive in V (coercive 
in a non reflexive Banach space), it may have and may have no solution, 
which means an appropriate regularization is needed. We do not investigate 
this issue further in a this work and only note that, after discretization in 
time, the problem becomes equivalent to a non-smooth optimization problem 
for each time layer. Indeed, let q = [u n+l — u n )/At and u = u n + ^q be 
approximate values at t — At(n + |). For each time layer we obtain 

m - <K<?) + (V • {u n + ^q}, V-{q-q})- {T n ^\ V • {q - q}) > 0, 
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which is equivalent to 

q«+l=arg min {^(V • q, V • q) + </>(q) + (V • u n - F n+1/2 , V • q)} 

q eV 

(27) 
Provided the latter problem has a solution, approximate value of u on the 
next time layer can be found as u n+1 = u n + Atq n+ z. 



2.2 Electric field in the Bean model for superconductors 



The h- and j-variational formulations of the Bean model can be used for ef- 
ficient computation of magnetic fields, current densities and, therefore, forces 
and their moments in various applications of type-II superconductors. How- 
ever, even if the current density is known, the electric field is not determined 
in a unique way by the assumption that the directions of e and j coincide and 
the current- voltage relation 

\j\<Jc \j\<Jc^ e = (28) 

or another constitutive law described by a monotone multivalued graph is 
satisfied. That is why calculating the electric field in a superconductor is gen- 
erally a non-trivial task. 3 In particular, this field is necessary to estimate 
the local AC loss e ■ j = j c \e\ that causes heating and thermal instability of 
superconductors. 

For some simple configurations, the electric field in superconductors has been 
considered in [31,32] and, recently, for a generalized Bean model in [33]. Here 
we propose a completely different approach based not on the determination of 
the magnetic field and subsequent integration of Faraday's law along the flux 
penetration streamlines but on a direct variational reformulation of the Bean 
model in terms of electric field. 

Let W = H(curl; Q) = {<p G L 2 (Q) \ V x cp e L 2 (Q)} be the space of electric 
fields in the superconductive domain Q C M 3 . Assuming directions of j and 
e coincide and the Bean model relations are satisfied in this domain, for any 
test field e e W we obtain 

(V x h) ■ (e-e) = j ■ (e - e) < j c \e\ -je = j c \e\ - j c \e\. 

Integrating over Q we get (V x h, e — e) < 0(e) — 0(e), where 0(e) = j c J n \e\. 



3 The case of an infinite cylinder in a perpendicular external magnetic field is an 
exception. 
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On the other hand, 



(V x h, e - e) = (h, V x { e - e}) + f(hx{e- e}) 



n. 



Let u = —Aq + Jq edt, where the vector potential Aq = G*(j+j e )\t=o satisfies 
/i = Vx A . Then 



d t u = e, V x u = —h + / V x edt = —h — / d t hdt = —h. 

Jo Jo 



V x edt = —ho — I 
/o Jo 

The following inequality is thus obtained: for any e G W 



(Vx«t,Vx{e- d t u}) - f(hx{e- d t u}) ■ n + 0(e) - <j)(d t u) > 0. (29) 

In the general case this is not yet the needed e-formulation because it contains 
the tangential component of magnetic field on T. However, for an infinite 
cylinder in a parallel uniform external magnetic field h\r = h e (t), where h e is 
the field generated by the external current, so 

(Vxii,Vx{e- d t u}) - h e (t) ■ f n x (e - d t u) + 0(e) - 4>{d t u) > (30) 

and it is not difficult to see that, after discretization in time, (30) becomes 
equivalent to a non-smooth optimization problem similar to (27). 

Let us now return to the general case and consider an auxiliary boundary 
value problem in the exterior domain u, 4 

d t h + V x e = 0, V x h = j e , 

(31) 

h\t=o = ho, e T \ r = £, 

where e T = n x (e x n)|r and £ is a tangential field given on T. We choose 
the vector potential of external current as A e = G * j e and define the field 

/i e = Vx (A e — A e \ t=0 ) + h . In uo this field satisfies 

V X h e = j e , V • h e = 0, h e \t=0 = h . 

Let us set h = h e + H, e = —d t A e + E and rewrite problem (31) as 

d t H + Vx E = 0, V x if = 0, 

(32) 

ff \t=o — 0, E T \ r — Ei, 



Since the displacement current is omitted, in our model the electric field in an 
insulator (the outer space) is not unique. The magnetic field, however, is. 
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where £\ = £ + d t A e r . Defining C/ = J * £^rft and integrating in time the 
equations from (32) containing E yield 

H + V xU = 0, VxE = 0, 

(33) 

U T \ r = U, 
where VI = J ' £dt + A £tT (x, t) — A £jT (x, 0). 

Let us note that we also have V • H = a.e. in a;, so it makes sense to seek a 
solution H in the Sobolev space H l {uo). Denote X = {?/> e H 1 ^) | V x t/> = 
0}. The problem (33) admits a weak formulation, 

Find HeX such that (if, i/>) = I ' (U x i/>) ■ n, V^ E X, (34) 

where the normal n is directed inside u and the integral on the right is under- 
stood as the duality pairing on H' 1 ^ 2 x H 1 ^ 2 . Existence of a unique solution 
to this problem follows from the Lax-Milgram theorem. This defines a linear 
operator K : U -»■ JEf r | r acting from if-^^r) to H l / 2 (T). It is not difficult to 
see that this operator is symmetric in the following sense: for any two functions 
v, w G H~ 1/2 (T) there holds 

I (K(v) x w) ■ n = f (K(w) xv)-n. (35) 

Indeed, let H v and H w be the solutions of (34) with U = v and 1A = w, 
correspondingly. Then 



[H V ,H W ] 



f(vx H w ) n= f(vx K(w)) ■ n, 
- (w x H v ) ■ n = (w x K(v)) ■ n 



and (35) is proved. We also see that J r (v x K(v)) ■ n = (H v , H v ) > for any 

veH- l / 2 (T). 

Let us now choose £ to be the continuous tangential component of the electric 
field on T, £ = e T . Then 

U = / e T dt + A er (x, t) - A eT (x, 0) = u T + A 0r + A eT - A eT \ t=0 
Jo 

= U T + G* (j\t=0+je) 

and the tangential component of the magnetic field can be presented as 

h T = h eT + H T = K(u T ) + T, 
where T = h eT + K({G * (j\ t=0 + j e )} T ). 
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We can now rewrite (29) as a mixed variational inequality: 

Find u : [0, T] — ► W such that for any e e W and almost all t 

(V x u, V x {g - d t u}) - / r (A"(ttr) x {g - d t u}) ■ n- 

(36) 
Jr(:F x {g - dtu}) ■ n + 0(g) - 0(9 t u) > 

and u\t=o = — ^4.0- 

Discretizing in time as above and making use of the equality (35) we obtain, 
for each time layer, a non-smooth optimization problem: 

e n+l/2 _ ar g m j n $( e ) 

eeW 

where 

$(e) = — (V x e, V x e) + — / (e r x A"(e r )) • n + 0(e) 
4 4 Jr 

- At | ({tf(e?) + ^ l+1 / 2 } x e r ) ■ n + (V x u n , V x e) 

and w n+1 = w n + Ate n+1 / 2 . 



Conclusion 



To derive quasistationary critical-state models of the spatially extended dissi- 
pative systems considered above, we needed to specify only the possible direc- 
tion of system's evolution, and which changes of external conditions make the 
state of the system unstable. The rate with which such a system driven by the 
external forces rearranges itself is determined implicitly by some conservation 
law coupled with a condition of equilibrium. This rate appears in the model 
as a Lagrange multiplier related to the equilibrium constraint. Although the 
conservation laws and conditions of equilibrium may vary, the multiplicity 
of metastable states, typical of many dissipative systems, is usually a conse- 
quence of a unilateral constraint. This makes variational inequalities a suitable 
tool for modeling these systems. 

Typically, some of the physically relevant conjugate variables are eliminated 
in transition to the variational formulation of critical-state problems. It is, 
however, possible to derive dual formulations (mixed variational inequalities) 
in terms of these variables. Although arising mathematical problems need 
further investigation, we believe these dual formulations will also serve a basis 
for efficient numerical simulations. 
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